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!  *  ■  I.  INTRODUCTION 

•  *■  /  • 

i  ,  •  i  . 

i 

The ' stability  of  thin  liquid  films  adjacent  to  slender  bodies  in 
high  spe$d  gas  flows  is  examined.  Maintenance  of  such  films  on  the 
external  surface  .of  hypersonic  vehicles  is  highly  desirable  from  the 
standj?oint,  of  thermal  protection  from  aerodynamic  heating  during  the 
re-entry  from  interplanetary  missions.  Such  a  thin  cooling  film  may  be 
established  by  allowing  liquid  injected  in. the  stagnation  region  of  the 
body  to.be  swept  downstream  covering.. the  surface.  The  coolant  mass  flow 
rate  required  to  effect  protection  of  the  vehicle  will  depend  on  the 

*  •  i  •  ' 

surface  area  and  on  mass  loss  due  to  droplet. entrainment  and  vaporization. 

,  i 

The  interaction  of  the  gas-liquid  interface  leads  to  the  formation  of 

isurface  waves  in  the  liquid  film.  Depending  on  the  flow  conditions, 

these  waves  will  either  be  damped  (stable  film)  or  will  grow  with  time 
i  i’  ’  i 

.(unstable  film) .  In  the  event  that  the  wave . amplitudes  are  not 

"  i  ■ 

sufficiently  decreased  with  time,  droplet  entrainment  will  become  'the 

»  * 

dominant  mechanism  of  mass  loss  to  an  extent  that  an  integral  layer  of 
Ijiquid  will  no  lohger  bp  maintained  on  the  surface. 

Couplediwith  the  effect'  on  the  liquid  film,  the  presence  of  waves 

i  1  ■ 

on  the  surface  also  affects  the  gaseous  boundary  layer.  The  flow  of  a 
ga3  over  a  wavy  or  “pebbled"  surface  is  similar  to  flow  <A?exr  a  rough 

surface.  These  conditions  give  rise  to  an  U,  the  transport  of 

\ 

mass,  momentum,  and  energy  in  the  boundary  layer  and  hence  an  increase 

I  » 

*  i 

in  heat  transfer  and  shear  stress.  The.  evaporative  mass  transfer  may  be 
!  • 

considered  to  consist  of  the  sum  of  that  contribution  which  would 

i 

result  i£  the  gas-liquid  interface  were  perfectly  smooth  and  the 
additionai  contribution  due. to  the  effect  of  the  surface  roughness. 

Gater  and  b'Scuyer*  have, indicated  that  this  latter  effect  is  independent 

1  '  i 

of  the  m$se  flow  rate.1  This,  however,  ignores  the  fact  that  the  surface 

i 

roughness,  wliicfe  is  related  to  the  wave  amplitude,  may  be  a  function  of 

it 

ii.  t 
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the  mass  flow  rate.  Thus,  the  effect  or  changes  in  the  mass  flow  on  the 
behavior  of  the  interface  cannot  be  disregarded  even  in  attempts  to 
predict  mass  transfer  due  to  evaporation.  . 

The  liquid  film  cooling  phenomenon  has  been  characterized  by  the 
<?■.  .isr^nce  of  a  region  in  which  heat  is  either  absorbed  by  the  liquid 
iiim,  or  used  to  effect  a  change  to  the  vapor  phase,  and  a  region  in  which  • 
cooling  is  accomplished  by  heat  transfer,  to  the  vapor  in  the  boundary 
layer  after  the  film  has  been  totally  vaporized.  The  effect  of  an 
instability  occurring  at  some  point  on  the  surface  of  the  vehicle  would 
be  to  prematurely  bring  about  the  reduced  cooling  effect  associated  with 
the  total  vaporization  region. 

If  an  analysis  of  the  stability  problem  is  to  be  applicable’ to  a 
given  flow  system,  the  probable  instability  mechanisms  of  the  system 
must  be  known  so  that  equations  and  boundary  conditions  consistent  with 
these  mechanisms  may  be  formulated.  The  instability  mechanisms  are  of 
two  basic  types:  those  arising  due  to  energy  transfer  from  the  basic 
flow  of  the  gas  or  liquid,  as  in  Tollmien-Schlichting  instabilities  and 
those  arising  due  to  energy  transfer  from  the  perturbation  of  the  normal 
and  shear  stresses  induced  by  the  wave  formation,  such  as  the  pressure 
perturbation  (Kelvin-Helmholtz)  and  shear  perturbation  instability 
mechanisms. 

Waves  capable  of  exhibiting  the  former  type  of  instability  were 

2 

possibly  observed  by  Marshall  and  Saric  in  flows  of  very  thin  films  at 
relatively  high  liquid  Reynolds  numbers,  in  this  case,  the  wave  speed 
was  less  than  the  interface  velocity  (slow  waves)  which  is  a  necessary 
condition  for  Tollmien-Schlichting  instability  according  to  Miles  . 

Most  waves  observed,  however,  have  not  been  in  this  regime. 

Craik'  classified  the  types  of  waves  observed  with  an  incompressible 
external  air  flow  which  occurred  at  different  values  of  the  film 
thickness.  He  found  that  instabilities  may  occur  for  very  thin  films 


2 


while  they  do  not  occisr  for  thicker  films  despite  the  relatively  high 
damping  effect  at  low  Reynolds  numbers.  He  concluded  that  the  dominant 
mechanism  of  .these  instabilities  was  the  action  of  the  shear  perturbation 
which  was  in  phase  with  the  wave  slope  and  tended  to  pile  up  liquid  at 
the  wave  crests.  • 

If  the  external  air  flow  is  compressible,  the  instability  mechanisms 
related  to  the  pressure  disturbance  are  of  different  form  depending  on 
whether  the  flow  is  subsonic  or  supersonic.  For  subsonic  flow,  the 
disturbance  pressure  is  180°  out  of  phase  with  the  surface  wave 
amplitude;  hence,  minimum  pressure  occurs  at  the  crests  and  maximum  in 
the  troughs,  thus  drawing  liquid  to  the  crests  of  the  waves.  However, 
for  supersonic  flow,  the  pressure  perturbation  is  in  phase  with  the  wave 
slope.  Chang  and  Russell  and  hacht3heim0’ '  have  demonstrated  in  this 
case  that  the  dominant  mode  of  energy  transfer  between  the  external  flow 
and  the  perturbations  of  the  liquid  flow  is  through  supersonic  wave  drag. 

Chang  and  Russell  performed  the  stability  analysis  for  a  gas-liquid 

interface  including  the  effects  of  the  perturbations  imparted  to  the  gas 

motion.  However,  the  liquid  layer  they  considered  was  infinitely  deep 

and  initially  quiescent,  Nachtsheim  considered  the  type  of  waves 

8 

observed  by  Larson  and  Mateer  in  the  melt  layer  of  a  re-entry  vehicle* 

He  carried  out  the  analysis  for  three-dimensional  disturbances  in  thin 
films  adjacent  to  supersonic  flows  and  showed  that  the  equations  may  be 
reduced  to  two-dimensional  form  if  the  wave  number  and  Mach  number  are 
appropriately  defined.  In  addition,  the  experiments  of  Marshall  and 
Saric  suggest  that  two-dimensional  waves  occur  at  intermediate  Reynolds 
numbers.  Nachtsheim' s  work,  however,  was  primarily  concerned  with  the 
highly  viscous  melt  layer;  and  he  concentrated  on  solutions  for  vanishingly 
small  Reynolds  numbers,  in  addition,  the  results  he  produced  were  not 
consistent  with  any  particular  external  flow  but  rather  were  obtained 
by  allowing  the  non-dimensional  parameters  in  the  governing  equations 
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to  vary  independently.: 

The  present  investigation  considers  thin  -films  at  moderate  film 
Reynolds  numbers  with  high  speed  external  flows'-  typical  of  "conditions 
arising-  for  liquid; films  in  r  in try  environments.  While  the  free  stream 
Mach  numbers  in  such ; environments  are  generally  hypersonic,  both' 
supersonic  and  subsonic  gas  flows  are  Used  to  provide  the  perturbation 
'pressure  condition  at  the  interface,  the  former  corresponding  to  a 
boundary  layer  much  thinner  than  the  disturbance  wave  amplitude  and- the 
latter  to  a  boundary  layer  which  is  thick  compared  to  the  disturbance 
amplitude. 

For  the  analysis,  instability  mechanisms  which  couple  the  disturbance 
motion  of  the  gas  to  the  disturbance  motion  of  the-  liquid  are  included, 
while  those  involving  energy  transfer  from  the  basic  flow  through  the 
action  of  shear  perturbations  at  the  interface  or  through  the  action 
of  body  forces  (Rayleigh-Taylor  instabilities) ,  are  specifically 
excluded.  The  analysis  of  Nachtsheini  is  generally  followed  except  that 
the  equation  of  motion  is  arrived  at  from  a  two-dimensional  standpoint 
and  consideration  of  the  vanishingly  small  Reynolds  number  regime  is 
excluded.  Analytic  solutions  are  also  developed.  Results  based  on  the 
numerical  and  analytic  solutions  to  the  equation  are  presented. 

II.  STABILITY  ANALYSIS 

*  i 

(A)  The  Basic  Flow 
Liquid  Film 

For  the  flow  of  a  thin  liquid  film  over  a  flat  solid  surface 
adjacent  to  a  supersonic  gas  stream  and  in  the  absence  of  perturbations, 
a  very  nearly  linear  velocity  profile  will  be  established  provided  the 
film  is  sufficiently  thin.  The  equilibration  of  shear  stress  at  the 

l 

interface  imposes  a  velocity  on  the  surface  of  the  liquid  film  which  is. 
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held:  constant  in  time  and  space  precluding  any  energy  transfer  from  the 
basic  flow  to  the  perturbations., 

A . two-dimensional  cartesian  coordinate  system  is  introduced  such 
that  the  x-axis  lies  in  the  plane  of  the  interface  and  the  y-axis  is 
perpendicular  to  that  plane  and  pointing  away  from  the  solid  wall  (see 
Fig.  1) .  The  basic  flow  of  the  liquid  is  then  given  by 

u-U^d  +  y/h) 

v  =  0  -h  <  y  <  0  (1) 

P  -  p©~  p^gy 

where  U  ,  the  surface  velocity,  is  subject  to  the  condition 

*v 

-  Tg!y=0  (2) 

The  gravity  vector  here  is  assumed  to  act  in  the  negative  y  direction. 
This  is  consistent  with  the  linear  from  of  the  velocity  profile.  It 
should  be  pointed  out,  however,  that  for  thin  films,  the  deviations  from 
the  linear  profile  will  be  small  when  the  gravity  vector  is  arbitrarily 
directed.  Hence,  the  neglect  of  the  error  associated  with  the 
application  of  the  following  analysis  to  the  general  case  is  consistent 
with  the  approximation  in  the  subsequent  linearization  of  the  governing 
equation.  However,  it  must  be  stipulated  that  the  gravity  vector  has 
no  component  in  the  positive  y  direction,  in  order  to  preclude 
instabilities  of  the  Rayleigh-Taylor  type. 

Externa 1  Gas  Flew 

W— >i^MI  im—  ■«]  ■■ 

The  basic  flaw  in  the  gas  is  uniform  and  given  by 

u  =  Ue 
v  =  0 

•  MM 

p  =  pe 

P  =  pe 

The  boundary  layer  is  assumed  to  be  very  thin  compared  to  the  liquid 


film  depth  and.  its  effect  is  confined  to  the  satisfaction  of  Eq.  (2) » 

The  value  of  the  shear  stress  at  the  interface  may  be  determined  from  a 
boundary  layer  analysis. 

(B )  Disturbance  Motion 

Liquid  Film 

A  component  of  a  general  harmonic  disturbance  is  now  introduced  to 
the  motion  of  the  liquid  film.  The  wave-like  form  of  this  disturbance 
is  illustrated  in  Pig.  1.  If  the  x-axis  is  directed  perpendicular  to  the 
wave  fronts,  perturbations  normal  to  the  x~y  plane  will  vanish  and  the 
total  velocity  and  pressure  are  given  by 

u  =  Ut(lry/h)  +  u(y)e1(kx~a,t) 
v  =  v(y)ei(,cx-"t)  (4) 

p  *  R.-pigy  +  p(y)ei(k!t~"lt> 

The  wave  number  k  is  assumed  to  be  a  real  dimensional  constant, 
while  the  parameter  uj  is  taken  to  be  complex.  The  wave  length  of  the 
disturbance  in  the  direction  of  the  wave  normal  is  2u/k.  The  real  part 
of  w  is  the  frequency  of  the  disturbance  and  the  imaginary  part  is  the 
time  amplification  factor.  Hence,  the  criterion  for  neutral  stability 
(waves  that  do  not  grow  or  decay  with  time)  is  that  the  imaginary  part 
of  w  vanish. 

External  Gas  Flow 

A  similar  disturbance  is  introduced  into  the  gas  flow 
a  =  oe  +  u(y)ei(kx-“tl 

7  =  v(y)ei<1'x-“*>  (5) 

P  =  pe  +  P(y)ei(kx-M:> 

7=  pe+p(y)ei(kx-”'t) 
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(C)  Disturbance  Equations 


Liquid  Film 


The  equations  governing  the  disturbance  amplitude  functions  are 
obtained  by  introducing  the  velocities  and  pressure  including  the 
perturbation  terms  into  the  Navier-Stokes  and  continuity  equations  and 
disregarding  non-linear  terms  in  the  perturbation  quantities  on  the 
assumption  that  the  disturbances  are  small. 

-ituu  +  U,  (1+y/h)  iku  +  vU  ./h  =  -~p  +  v .  (u/,-k2u)  (6) 

-o  "t 

-iujv  +  U  .  (1+y/h)  ikv  =  — '  +  v.  (v"  -k3v)  (7) 

iku  4  v '  =  0  (8) 

Primes  here  indicate  differentiation  with  respect  to  y. 

Eliminating  p  then  u  from  Eq.  (7)  leads  to  the  Orr-Sommerfeld 
equation  for  the  liquid  film 

vlv-2k3v"  +  k4v=  ~|[U^ (1+y/h) -wAJ  (v"~k3v) } .  (9) 

and  eliminating  u  from  Eq.  (6)  leads  to  the  following  expression  for  the 
disturbance  pressure  in  terms  of  v  and  its  derivatives 

P  =  -^{iC^-^U+yANJv'+v  -^3+  ^(v"-k3v')|  (10) 

The  disturbance  shear  stress  in  the  film  is 

Txy  =  kv"+  ikv)  (11) 

E:iternal  Gas  Flow 

The  equations  of  motion  for  the  gas  flow  are  derived  in  a  similar 
manner.  The  velocities,  pressure  and  density  are  introduced  into  the 
inviscid  momentum  equations  and  the  continuity  equation  for  a  compressible 
fluid.  Non-linear  terms  are  again  neglected. 


(Ue-u>A)u  =  ~ae  P 
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(13) 
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P  ,  f -^vT-S' vTJSIK^JSSfT*; 


and  the  kinematic  condition  is 

,  V  =  -^  +  u 

v  at  u  ax 

In  terms  of  perturbation  quantities 


v  =  ikU  e 
e 


The  solution  to  Eq.  (17)  subject  to  boundary  condition  (20)  is 

v  =  ikUe«eifcAea-1  y  (21) 

Now  c  may  be  determined  from  the  kinematic  conditiai  for  the  liquid 


at  the  interface 


v  =  ik(U^-uJA)c 


Then  the  disturbance  pressure  in  the  gas,  Eq.  (16)  reduces  to 

PeoVv 


P  =  (U.-ojA)^' ^“1 


in  the  supersonic  case  and  to 


ip.Vv 

(u^-u)A)vT:ir3' 


in  the  subsonic  case. 


Thus,  the  condition  on  the  normal  stress  at  the  interface  in  the 
supersonic  and  subsonic  cases  are  given,  respectively*  by 


i  v.  U. 

vw+  [{U^-u)A)-3iv^k3v' --jp  v 

.  ikp  U  s 

«  f-(Tk3+p  ,g+~d =)  (  Vk  ) 

H  *  Vm_s-i  u i  /K 


iv  U 

vw+  C(U^-«JA)-3iv^k]v ' — v 

i  ^  Pe^e 3  v 

■  -±-(Tk»+p  g-  77T~w 

H  1  JT=m^  w/ 


If  the  disturbance  shear  stress  in  the  gas  is  neglected,  the 
boundary  condition  becomes 
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I 


v*  +  k8v  =  0 


at  y  =  0 


The  condition  at  the  wall  is-that  the  velocity  vanish 


v  -  v'  =  0 


at  y.  =  -h 


(E)  Dimensionless-  Form  of  Equations 


The  variables  may  be  non-dimensionalized  as  follows 


r\  =  y  /h 

*  (n  )  *  v(y)/U^ 

Equation  (9)  and  the  boundary  conditions  (24) ,  (25) , (26)  and  (27) 


are  transformed  to 


at  t]  =  0 


at  r\  =  -1 


$iv-2aa*"+a4*  =  iaRCd+n-c)  (*"~as*)  ] 


$  "  +  a2  $  =0 


§  =  *'  =  0 


where  primes  now  denote  differentiation  with  respect  to  r)  and 

g/k 

c  ss 


a  =  kh 


V  “ 


c  - 

f“  hpeV 

F  =  i 


H  =  -  (supersonic) 

cfRvMe^r 


CfiVl^7 


(subsonic) 
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If  the  function  cp  is  defined  such  that 

p(Ti)  =  §(*)/§  (0)  =  v(y)/v(0)  (32) 

the  equations  and  boundary  conditions  will  retain  their  same  form  as 
follows: 

cp^ 2 as ep /y  +  a4cp  =  iaR[(l+TT~c)  (cp"-asco)  }  (33) 

at  n  =  0 

•^co^+Cd-c)-  “~]cp'~cp  =  (gr  +  p-+H)  ./  (34) 


<$"  +  aacp  =  0  _  ...  (35) 

at  r)  =  -1 

<p  =  cp'  =  0  .  (36) 

with  the  additional  condition  that 

cp  *  1  at  n  =  0  (37) 

Noting  that  the  combination  cp^  -antp^n""2'  appears  repeatedly  in  the 
equation,  a  new  function  iji  may  be  defined 

^  -  a3co  (38) 

The  equation  may  now  be  reduced  to  a  system  of  two  second-order 

•equations 

cp"  =  aBcp  + 

(39) 

i|ftf  -  a3i)  =  iaR(l+iT-c)  i) 

and  the  boundary  conditions 
at  t]  =  0 

t'  =  2a3cp/-ia{[R(|J+  ^-+H)3  -^.-Rttc-Dcp'+l]}  (40) 


and  at  r)  =  -1 


(f  +  2ascp  =  0 

cp  =  cp'  =  0 


(41) 

(42) 
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The  above  equations  and  boundary  conditions  constitute  an  eigenvalue 
problem.  .The  eigenvalue  consists  of  Refc},  the  non-dimensionalized  wave 
speed  and  Imjcj, ‘the  time  amplification  factor.  If  the  eigenvalue  is 
correctly  chosen,  the  homogeneous  boundary  condition  (40)  will  be 
satisfied  independent  of  the  choice  of  the  initial  value  of  cp.  This  is 
illustrated  by  the  fact  that  the  equations  retain  the  same  form  under 
the  transformation  (32)  which  is  equivalent  to  the  choice  of  initial 
condition  (37) . 

For  a  given  liquid  and  gas,  the  problem  is  completely  specified  by 
the  geometry,  mass  flow  rate,  and  the  balance  of  shear  stress  for  the 
basic  flow.  Geometric  considerations,  in  general,  lead  to  an  expression 
for  the  mass  flow  of  the  form 


(F;  Analytic  Solution  for  Small  Values  of  g 


If  the  boundary  condition  (40)  is  to  remain  well  behaved  at  c=l, 
then  a  must  be  equal  to  zero  for  c-1.  Hence,  solutions  are  sought  for 
vanishingly  small  a  with  c  near  1.  The  parameter  c  and  the  dependent 
variables  cp  and  are  expressed  as  a  series  in  powers  of  a 


c  =  1  +  c:  a  +  csa3  +  c3  a3  +  . 
cp  =  cpQ  +  cpi  a  +  cps  as  +  cp3  a3  + 
=  »  +  Vi  a  +  a8  +  fs  a3  + 


Substitution  into  the  differential  equations  yields  to  zeroth 
order  in  a 
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with  the  boundary  conditions 
at  n  =  0 


at  t)  =  -1 


=  1 
*o  "  0 


i  /  iR 

^o  cxP* 


<P0  =  0 


•o'  =  0 


the  solution  for  co  and  *  are 

o  To 


<»0  =  “3Sti3+3/s  T)  +  1 


*0  =  -3r, 


and  from  the  extra  condition  at  r|  =  0 

„  _  -iR 

ci  =  3pr 


(44) 


(45) 


(46) 


(47) 


retaining  first  order  terras  in  o  leads  to  the  following  equations  and 
boundary  conditions  in  the  supersonic  case 


If  the  same  procedure,  is  carried  out  in  the  subsonic  case,  the  resulting 


expression  for  c  is 


.  .  iR  «  ^  /  7R?  , 

l-3pTU+  (72PS  + 


acfViiMp- 


:).a8 


(58) 


,  '  ^  r  /  1  .  1  .  2  ,R3  1  ,  13  ■,  ■  7R  , 

+  lR[(4536  +  36  +  45Fy)F3'  '  3w5"f30F5'^  ~72c£VT^Meg'  a 

with  no  change  in  the  expression  for  cp. 

»  I 

Note  that  the  expression  foi:  c  indicates  that  a=0  corresponds  to 
neutral  stability  in  either  case.  - 


(G)  Exact  Subsonic  Solution 

,  In  the  subsonic  case  a  simple  analytic  solution  giving  the  cut-off 
wave  numbers  (for  which  lm{c|=0)  may  ne  obtained  when  neutral  stability 
implies' c=l.  In  the  subsonic  case,H  is  real  and  negative.  Hence,  the 
term  in  brackets  multiplying  l/(c~l)  in  the  boundary  condition  (40)  may 

;  I  i 

vanish  and  must  do  so  as  c  approaches,  unity.  Therefore,  the  cut-off 

1 

I 

wave  nuijxbers  may  be  determined  from  the  expression  ' 

2  '  . 


StL. 

wa 


+  cTS*  “ 


°N 


:£R/L-M{ 
1 

The  following  values  of  are  obtained 
1  ' . 

■  .  1  V* 


=  0 


(59) 


W 


2cfR*/l-Mes 


J(. 

^l2cfRvT=M “ 


”  z 

\-ih) 


(60) 


(H)  Remarks  on  the  Applicability  of  Subsonic  and  Supersonic  Solutions  - 
Boundary  Layer  Effects  ~  *  '  '  “ 

While  the  present  investigation  is  concerned  exclusively  v;ith 
1  1  . 

situations  in  which  th6  external  gas  stream  is  supersonic,  the  subsonic 

solutions  may  also  be  of  interest.  Nachtsheim  has  suggested  that  the 

basic  flow  in  the  gas  may  be  regarded  as  a  "suitable  mean"  of  an  actual 

boundary  layeri  flow.  However,  the  presence  of  the  gaseous  boundary 

layer  will  influence  the  -phase  lag  of  the  pressure  perturbation  at  the 
1  ■  ' 
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1 


l 


l 


1 


interface  as  well  as  the  appropriate  values  of  the  Mach  number  and 
"free  stream"  velocity.  Indeed,  the  subsonic  solution  may  be  valid  for 
a  flow  with-  a  supersonic  external  stream  when-  the  boundary  layer  is 
thick  enough  to  attenuate  the  pressure  perturbations  in  the  subsonic 
region. 

In  the  limit  of  an  infinitely  thin  boundary  layer  with  a  supersonic 

external  stream,  the  pressure  perturbation  is  seen  to  exhibit  a  90°  phase 

lag  with  respect  to  the  surface  wave  amplitude  while  the  Mach  number 

appropriately  assumes  the  free  stream  value.  In  the  opposing  limit, 

as  Me  boundary  layer  thickness  becomes  large  compared  to  the  wave 

amplitude,  the  external  flow  may  be  regarded  as  subsonic  for  which  a 
o 

180  phase  lag  occurs.  Hence,  as  the  ratio  of  boundary  layer  thickness 
to  wave  amplitude  varies  between  the  two  limits,  it  is  to  be  expected 
that  the  phase  lag  will  vary  between  90°  and  180° •  This  behavior  has 

g 

been  amply  demonstrated  experimentally  by  Muhlstein  and  Beranek  . 

In  order  to  take  this  effect  into  account,  the  pressure 
coefficient  at  the  interface,  H,  could  be  expressed  as  follows: 


cfRv'|Mea-i| 


where  3  is  the  phase  lag.  This  clearly  reduces  to  the  supersonic  case 
when  0=90°  and  to  the  subsonic  case  when  0=180°.  Here  0,  cf  and  Mg 

’ff 

would  be  functions  of  6  /s.  The  difficulty  in  applying  this  approach 

lies  in  making  an  appropriate  choice  for  the  parctmeter  values.  The 

results  of  Muhlstein  and  Beranek  indicate  that  it.  is  sometimes  reasonable 

to  use  the  limiting  values  of  0  (particularly  in  the  subsonic  case)  when 
•/( 

6  /e,  cf  and  M  are  not  very  near  their  limiting  values. 
a*  y 

Therefore,  the  approach  here  will  be  to  employ  the  strictly 
supersonic  and  subsonic  phase  lags  in  the  analysis.  The  problem  of 
choosing  parameter  values  is  not  clearly  resolved,  although  reasonable 
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choices  may  be  made . 


IIX>.  -NUMERICAL  PROCEDURE 


S  * 

Since  the  parameters  R,  cf,  F,  W,  and  Me  are  specified  for  a 
given  liquid  mass  flow  and  given  external  flow#  it  is  necessary  to 

r  *  ‘  € 

determine  the  eigenvalue  c  as  a  function  of  a  satisfying  the  equations 
and  boundary  conditions.  Thus,  for  any  a,  a  numerical  procedure  may  be 
used  to  determine  the  values  of  c  for  which  Eq.  (39)  has  solutions 
satisfying- the  boundary  conditions  (40)-(42). 

In  order  to  carry  out  a  numerical  integration  of  the  equations 
starting  at  17=0,  an  initial  guess  for  the  eigenvalue  c  and  two  additional 
conditions  are  required.  The  initial  conditions  may  be  expressed  in  the 
form 


cp(0)  =  A 
cp'(O)  »B 


(61) 


As  previously  noted,  solutions  satisfying  the  boundary  conditions 
are  completely  independent'  of  the  choi.ce  of  the  initial  condition  on  cp. 
Therefore,  A  is  taken  to  be  1  as  in  the  analytic  solution.  The  value 
of  B,  therefore,  must  be  determined  from  the  boundary  conditions.  For 
small  values  of  a,  Eq.  (57)  differentiated  with  respect  to  17  and 
evaluated  at  17=0  may  be  used  to  obtain  an  initial  guess  for  B  and 
Eq.  (56)  or  Eq.  (58)  may  be  used  to  find  c.  The  formulated  initial 
conditions  (40)  and  (41)  may  then  be  written  in  terms  of  the  assumed 
initial  condition  as  follows: 


$  =  -?.a2 


=  2a*B  +  ia{[R(~-+p-4  H)  ]  —■  4  R[(c-1)B  4  13}  (62) 

The  equations  are  then  integrated  from  r\  -  0  to  17  -  “1  and  the 
solutions  are  checked  to  see  if  the  conditions  at  the  solid  surface  are 
satisfied.  Since  cp  and  cp'  must  both  vanish  at  n  =  ~l»  the  sum  of  their 


absolute  values*  constitutes  the  error  corresponding  to  any  guess  for  c 
and  B.  I-f  the  boundary  conditions  are  not  satisfied  to  within  the 
accuracy  of  the  integration,  c  and  B  must  be  adjusted  accordingly. 

In  order  to  facilitate  the  rapid  convergence*  to  a  satisfactory 
solution,  it  is  helpful  to  know  how  cp  and  cp'  vary  with  B  and  c  at  t|  =  ~l. 
Differential  equations  for  c^,  cpg ',  cpc  and  cpc ',  where  subscripts  indicate 
differentiation,  may  be  obtained  by  differentiating  the  Eqs.  (39)  and 

initial  conditions  (61)  and  (62)  with  respect  to  B  and  c.  The  resulting 

*« 

equations  and  initial  conditions  are: 


cPb"=o3«Pb+  *B 
i)b"  =  ria!)B  +  i nR ( 1+  T| -  c)  i)B 

Cpc/'  =  u^c+  *c 

>)c"  =  a3 +  icR(l  +  ri -  c)  i)c-ioRi|f  (63) 

%  “  0 

V  =  1 
"  0 

'  =  2on  -  iaR(c-:l)  (64) 

D 


Vc'“  0 


V  =  0 
*'c  =  0  .. 


‘  !)c,  =  -icx|C?(^r+^+H)]^rp+ SB}  (65) 

The  adjustments  AB  and  Ac  to  the  ith  guess  are  obtained  from  the 
Newton-Raphson  equations:  ■' 

cpi+i  (-1)  -  «Pi  (-1)  +  ^  .(-D  AB  +  cp^-l)  AC  =0 

cpi+i'  (-1)  =>^'(-1)  +  .CPg  7  (~1)  AB  +  cpc  '(-l)  Ac  =0 
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(66) 


A  program  incorporating  the  above  numerical  approach  was  run  on 
the  IBM  360/50  computer  at  >the  Polytechnic  Institute  of  Brooklyn  and 
on  -the -remote  terminal  of' the  CDC  6600  computer  located  at  the  Preston 

Jj 

Ri  Bassett  Research  Laboratory. 

The  numerical ^integration  was  carried  out  using  the  fourth  order 
Runge-Kutta  subroutine,  RKGS,  taken  from  the  IBM  system/360  Scientific 
Subroutine  Package.  This  subroutine  allows  the  programmer  to  specify 
the  tolerable  integration  error  and  adjusts  the  step  size  until  this 
specification  is  met.  In  order  to  make  use  of  this  subroutine it  was 
necessary  to  recast  the  equation  in  terms  of  real  variables  and  reduce 
them  to  first  order.  This  led  to  a  system  of  twenty-four  first  order 
linear  differential  equations.  The  eigenvalue  and  initial  conditions 
were  adjusted  by  the  Newton-Raphson  technique  until  the  error  defined 
by  the  sum  of  absolute  values  of  cp(-l)  and  cp'(-l)  was  zero  to  within  the 

_E 

accuracy  of  integration.  This  usually  meant  less  than  5  x  10  . 

After  this  technique  had  converged  to  an  eigenvalue,  a  was 
incremented  and  the  corresponding  new  eigenvalue  was  obtained.  This 
procedure  was  continued  using  a  second  degree  extrapolation  to  predict 
the  new  eigenvalue  and  initial  condition. 

It  was  found  that  an  increasingly  accurate  initial  guess  was 
required  to  insure  convergence  to  the  proper  eigenvalue  as  a  was 
increased.  As  a  consequence,  it  was  necessary  to  keep  the  increment  of 
a  very  small.  A  highly  accurate  initial  guess  is  also  desired  since  the 
integration  procedure  is  time  consuming  and  less  rapid  advancement  in 
a  is  preferable  to  permitting  numerous  iterations  of  the  integration. 

In  addition,  the  possibility  of  convergence  to  an  extraneous  eigenvalue 
exists  if  care  is  not  taken  in  selecting  the  first  approximation  to  the 
eigenvalues  and  initial  conditions.  The  behavior  of  the  eigenfunctions 
for  two  values  of  a  is  illustrated  in  Fig.  2.  Typically,  the  curvature 
of  these  functions  increases  with  increasing  a  and  thus  the  computation 
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time;  become a  .significantly/  greater,  in  some  cases  this  curvature  was 
sufficiently  large  .that  an  inordinately  small  step  size  was  required  and 
the  integration  procedure  failed.  When  these  cases  were  rerun  with  a 
relaxation  in  the  accuracy  criterion,  such  that  the  solution  could  be 
.obtained,  it  appeared  that  the  .eigenvalues  obtained  were  relatively 
insensitive  to  the  accuracy  requirement. 

IV.  PRESENTATION  AND  DISCUSSION  OP  RESULTS 
(A)  Flow  Conditions 

The  external  flow  conditions  used  .are  given  in  the  following  table. 


mm 

H9B 

<WIV  -1 

2.00 

Supersonic  Case 

2.00 

.0115 

.020 

Subsonic  case 

0.25 

.0264 

The  Weber  and  Froude  numbers  consistent  with  a  surface  shear  stress 
of  .005  psi  are  given  by 

W  =  .04958R3/4 
F  =  23.793R1/4 

(B)  Supersonic  Case 

The  behavior  of  c  as  a  function  of  a  is  plotted  in  Fig.  3  for 
several  values  of  the  Reynolds  number.  There  are  characteristically 
two  cut-off  wave  numbers  other  than  a  =  0.  It  is  possible  to  use 
Eq.  (56)  to  determine  the  first  cut-off  wave  number,  ,  when  the 
coefficient  of  a3  is  positive. 

In  Fig.  4  the  cut-off  wave  numbers  are  plotted  versus  Reynolds 
number.  Most  striking  about  their  behavior  is  that  there  exists  a 
minimum  Reynolds  number  below  which  no  instabilities  occur  for  any  wave 
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number.  This  can  be  attributed  to  that  fact  that  for  small  Reynolds 
number  the  stabilizing,  effect  of  the  svurface  tension  becomes  dominant. 
That  is,  in  the  boundary  condition  (40)  the  term  a2/M2  dominates  for 
small  n.  Further,  since  the  stabilizing  term  1/FS  dominates  for  large 
R,  it  might  be  expected  that  there  exists  a  maximum  Reynolds  number 
above  which  all  wave  numbers  are  stable.  However,  it  has  not  been 
possible  to  demonstrate  that  here.  Results  obtained  in  previous 
publications  which  ignored  the  necessary  R  dependence  of  F  and  W  have 

i  - 

predicted  such  a  maximum  Reynolds  number,  but  these  results  must  be 
regarded  as  inapplicable  to  real  systems.  Flows  for  which  there  are  no 
values  of  a  for  which  lm{c|>0may  be  regarded'  as  generally  stable. 

(C)  Subsonic  Case 

Comparison  of  the  analytic  and  numerical  solutions  indicates  that 
the  former  is  exact  for  all  values  of  R.  The  numerically  obtained  plot 
of  amplification  factor  versus  wave  number  is  displayed  in  Fig.  5  from 
which  the  two  characteristic  cut-off  wave  numbers  are  evident.  In 
addition,  the  wave  speed  is  seen  to  be  smaller  than  unity  in  the  unstable 
region.  This  behavior  is  also  clear  in  Fig.  6  where  c  is  plotted  in 
the  complex  plane.  For  all  values  of  R,  c=l  at  neutrally  stable  values 
of  a,  c  <  1  at  unstable  values  of  a  and  c>l  at  stable  values  of  a.  hs 
a  is  increased,  starting  at  the  first  cut-off,  the  branches  are  traversed 
in  a  counter-clockwise  manner.  This  effect  rather  than  the  Tollmien- 
Schlicting  instability  mechanism  may  explain  the  slow  waves  observed  by 
Marshall  and  Saric.  It  is  notable  that  these  waves  were  observed  only 
in  the  thinnest  films  considered  for  which  the  ratio  of  boundary  layer 
displacement  thickness  to  wave  amplitude  was  greatest. 

The  a-R  diagram  presented  in  Fig.  7  indicates  that  no  finite 
maximum  or  non-zero  minimum  Reynolds  numbers  exist  beyond  which  there 
are  no  unstable  wave  numbers.  This  is  generally  true  for  values  of 
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cf  below  some  maximum  as  indicated  in  Fig.  8.  Thus,  such  flows  are 
generally  unstable  when  c^  lies  below  the  maximum  and  stable  when 
lies  above  the  maximum. 

V.  CONCLUSIONS 

The  following  general  conclusions  may  be  drawn  from  the  present 

results  when  the  pressure  perturbation  instability  mechanism  is 

predominant.  For  some  values  of  the  Reynolds  number,  no  instabilities 

will  occur  when  the  phase  lag  is  180°.  For  boundary  layer  flows  where 

o  ° 

0  is  less  than  180  the  effect  of  approaching  0=90  may  be  regarded  as 
stabilizing.  As  the  unstable  slow  waves  grow  in  the  boundary  layer,  the 
effect  is  to  retard  the  appropriate  phase  lag  and  stabilize  the  motion 
for  certain  values  of  R.  In  this  manner,  steady  amplitude,  fast  or  slow 
wayes  may  be  produced  without  non-linear  effects. 


22 


VI .  REFERENCES 

1.  Gater,  R.A.  and  L'Ecuyer,  M.R.s  A  Fundamental  Investigation  of  the 
Phenomena  that  Characterize  Liquid  Film  Cooling,  Purdue  University 
Report  TM-69-1,  January  1969. 

2.  Marshall,  B.W.  and  Saric,  W.S.:  An  Experimental  Investigation  of 
the  Stability  of  a  Thin  Liquid  Layer  Adjacent  to  a  Supersonic 
Stream.  AIAA  Paper  No.  70-801,  1970. 

3.  Miles,  J.W.s  The  Hydrodynamic  Stability  of  a  Thin  Film  of  Liquid 
in  Uniform  Shearing  Motion.  J.  Fluid  Mech.,  Vol.  8,  Part  4, 

pp.  593-610,  1960. 

4.  Craik,  A.D.D.:  Wind-Generated  Waves  in  Thin  Liquid  Films.  J. 

Fluid  Mech.,  Vol.  26,  Part  2,  pp.  369-392,  1966. 

5.  Chang,  I.D.  and  Russell,  P.E.:  Stability  of  a  Liquid  Layer 
Adjacent  to  a  High  Speed  Gas  Stream.  Phys.  Fluids,  Vol.  8,  No.  6, 
pp.  1018-1026,  1965. 

6.  Nachtsheim,  P.R.:  Analysis  of  the  Stability  of  a  Thin  Liquid  Film 
Adjacent  to  a  High-Speed  Gas  Stream.  NASA  TN  D-4976,  January  1969. 

7.  Nachtsheim,  P.R.s  Stability  of  Crosshatched  Wave  Patterns  in 
Thin  Liquid  Films  Adjacent  to  Supersonic  Streams.  Phys.  Fluids, 
Vol.  13,  No.  10,  pp.  2432-2447,  1970. 

8.  Larson,  H.K.  and  Mateer,  G.G.:  Crosshatching  -  A  Coupling  of  Gas 
Dynamics  with  the  Ablation  Process.  AIAA  Paper  No.  68-670,  1968. 

9.  Muhlstein,  L.  and  Beranek,  R.G.:  Experimental  Investigation  of  the 
Influence  of  the  Turbulent  Boundary  Laver  on  the  Pressure 
Distribution  Over  a  Rigid  Two-Dimensional  Wavy  Wall.  NASA  TN 
P-6477,  August  1971. 

10.  Betchov,  R.  and  Criminale,  W.O.:  Stability  of  Parallel  Flows. 
Academic  Press,  New  York,  1967. 

11.  Liepmann,  H.W.  and  Roshko,  A.:  Elements  of  Gas  Dynamics.  John 
Wiley  and  Sons,  Inc.  New  York,  pp,  208-215,  1967. 


23 


FIG.  (I )  FLOW  CONFIGURATION 
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BEHAVIOR  OF  THE  EIGENFUNCTION'. 
SUPERSONIC  CASE 
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FIG.  (2b)  BEHAVIOR  OF  THE  EIGENFUNCTION, 
SUPERSONIC  CASE 
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FIG. (7)  GUT -OFF  WAVE  NUMBERS  VS.  REYNOLDS 
NUMBER,  SUBSONIC  CASE 


